﻿
using System;
using System.Threading;

namespace LxCFD
{
    public partial class TurbulenceModelkw : ScalarHeadBase
    {
        public static readonly double Clim = 7 / 8.0;
        public static readonly double alpha = 13 / 25.0;
        public static readonly double betastar = 9 / 100.0;
        public static readonly double div_2_betastar = 2 / betastar;
        public static readonly double sigma = 1 / 2.0;
        public static readonly double sigmastar = 3 / 5.0;
        public static readonly double sigmado = 1 / 8.0;
        public static readonly double betao = 0.0708;
        public static readonly double Cmu14 = Math.Sqrt(Math.Sqrt(0.09));
        public static readonly double Karman = 0.41;
        public static readonly double Karman_Reciprocal = 1.0 / Karman;// 1 / 0.41;
        public static readonly double B = 5.25;
        public static readonly double beta1 = 0.075;


        public NSEquation k;
        protected double[] k_x;
        protected double[] k_b;
        protected double[] k_A_Diag;
        protected Vector3d[] k_grad;
        protected double k_minLim;
        public NSEquation w;
        protected double[] w_x;
        protected double[] w_b;
        protected double[] w_A_Diag;
        protected Vector3d[] w_grad;
        protected double w_minLim;
        protected double[] w_A_LU_O;
        protected double[] w_A_LU_N;


        public Value1FromArray1<double> Gamma;
        protected double[] GammaCommon;
        protected double[] cellCommon;//边界处储存y+值，内部紧邻边界处存储Pk或omegaC值，内部存储临时的k方程或w方程Gamma值

        protected VelocityVector VelVector;
        protected Vector3d[] U_grad;
        protected Vector3d[] V_grad;
        protected Vector3d[] W_grad;
        protected Value1Base<double> density;
        protected Value1Base<double> dynVis;

        protected double[] WallDis_x;
        protected bool[] isWallAdj;


        public override void CheckCreate(FvZone zone)
        {
            base.CheckCreate(zone);
            if (k == null)
            {
                k = new NSEquation();
                k.CheckCreate(zone);
            }
            this.k_x = k.x;
            this.k_b = k.b;
            this.k_A_Diag = k.ADiag;
            this.k_grad = k.Grad;
            if (w == null)
            {
                w = new NSEquation();
                w.CheckCreate(zone);
            }
            this.w_x = w.x;
            this.w_b = w.b;
            this.w_A_Diag = w.ADiag;
            this.w_grad = w.Grad;
            this.w_A_LU_O = w.ALUO;
            this.w_A_LU_N = w.ALUN;
            if (cellCommon == null)
            {
                cellCommon = new double[zone.Cells.Length];
            }
        }
        public override void Initialize()
        {
            k.Initialize();
            // MultiTask.Sync(taskIndex);
            k.Limit_x();
            // MultiTask.Sync(taskIndex);
            w.Initialize();
            //MultiTask.Sync(taskIndex);
            w.Limit_x();
        }
        public void SetZone_part1(VelocityVector velocityVector, Value1FromArray1<double> gamma, Value1Base<double> density, WallDistance WallDis)
        {
            this.Gamma = gamma;
            this.GammaCommon = gamma.CommonArray;

            this.VelVector = velocityVector;
            this.U_grad = velocityVector.U.Grad;
            this.V_grad = velocityVector.V.Grad;
            this.W_grad = velocityVector.W.Grad;

            this.density = density;
            this.dynVis = velocityVector.DynamicViscosity;

            this.WallDis_x = WallDis.x;
            this.isWallAdj = WallDis.isWallAdj;
        }
        public void SetZone_part2(ValueBase<double> k_S, Value1Base<double> k_iniValue, double k_relax, ConvectionSchemes k_convect, double k_minLim, ValueBase<double> w_S, Value1Base<double> w_iniValue, double w_relax, ConvectionSchemes w_convect, double w_minLim)
        {
            k.SetZone_CommonTerm(Gamma, k_S, k_iniValue, k_relax);
            k.SetZone_ConvTerm(new Value1Const.Double(1), k_convect, VelVector);
            this.k_minLim = k_minLim;
            k.SetLimitRule(true, true, k_minLim, false, double.MaxValue);
            w.SetZone_CommonTerm(Gamma, w_S, w_iniValue, w_relax);
            w.SetZone_ConvTerm(new Value1Const.Double(1), w_convect, VelVector);
            this.w_minLim = w_minLim;
            w.SetLimitRule(true, true, w_minLim, false, double.MaxValue);
        }
        public void SetZone_part2(ValueBase<double> k_S, double TurbulenceIntensity, double k_relax, ConvectionSchemes k_convect, double k_minLim, ValueBase<double> w_S, double VisRate, double w_relax, ConvectionSchemes w_convect, double w_minLim)
        {
            BoundkValue boundk = new BoundkValue();
            boundk.SetParameter(TurbulenceIntensity, VelVector);
            k.SetZone_CommonTerm(Gamma, k_S, boundk, k_relax);
            k.SetZone_ConvTerm(new Value1Const.Double(1), k_convect, VelVector);
            this.k_minLim = k_minLim;
            k.SetLimitRule(true, true, k_minLim, false, double.MaxValue);
            BoundwValue boundw = new BoundwValue();
            boundw.SetParameter(VisRate, k_x, Cells, density, dynVis);
            w.SetZone_CommonTerm(Gamma, w_S, boundw, w_relax);
            k.SetZone_ConvTerm(new Value1Const.Double(1), w_convect, VelVector);
            this.w_minLim = w_minLim;
            w.SetLimitRule(true, true, w_minLim, false, double.MaxValue);
        }

        public void SetBound_WallOrSymmetry(FvThread thread)
        {
            k.SetBoundary_symmetry(thread);
            w.SetBoundary_symmetry(thread);
        }
        public void SetBound_Flow(double TurbulenceIntensity, double TurbRadio, FvThread thread)
        {
            BoundkValue boundk = new BoundkValue();
            boundk.SetParameter(TurbulenceIntensity, VelVector);
            BoundwValue boundw = new BoundwValue();
            boundw.SetParameter(TurbRadio, k_x, Cells, density, dynVis);
            SetBound_Flow(boundk, boundw, thread);
        }
        public void SetBound_Flow(Value1Base<double> k_flowValue, Value1Base<double> w_flowValue, FvThread thread)
        {
            k.SetBoundary_flow(k_flowValue, thread);
            w.SetBoundary_flow(w_flowValue, thread);
        }
        public void Cal_MomGamma(int taskIndex)
        {
            cal_muT_innerCell(taskIndex);
            MultiThreads.Sync(taskIndex);
            calGamma_byAverage(InnerIndex, taskIndex);
            foreach (var b in VelVector.NoSlipWalls)
            {
                b.ThreadIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
                for (; i1 < endp1; i1++)
                {
                    GammaCommon[i1] = calBound_mu_w(i1, ref Faces[i1]);//calBound_muw(i1);
                }
            }
            foreach (var b in VelVector.Symmetries)
            {
                calGamma_byBoundOwner(b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelVector.VelocityInlets)
            {
                calGamma_byBoundOwner(b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelVector.PressureInOuts)
            {
                calGamma_byBoundOwner(b.ThreadIndex, taskIndex);
            }
            //测试用，先直接选用层流参数
            //totalIndex.FaceIdx.GetSeg(out int i1, out int endp1, taskIndex);
            //for (; i1 < endp1; i1++)
            //{
            //    GammaCommon[i1] = dynVis.getFaceValue(i1);
            //}
        }

        public void Cal_kGamma(int taskIndex)
        {
            cal_kwVis(TotalIndex, sigmastar, taskIndex);
            MultiThreads.Sync(taskIndex);
            calGamma_byAverage(InnerIndex, taskIndex);
            foreach (var b in VelVector.VelocityInlets)
            {
                calGamma_byBoundNeighbor(b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelVector.PressureInOuts)
            {
                calGamma_byBoundNeighbor(b.ThreadIndex, taskIndex);
            }
            //测试用，先直接选用层流参数
            //totalIndex.getFaceSection(out int i1, out int endp1, taskIndex);
            //for (; i1 < endp1; i1++)
            //{
            //    GammaCommon[i1] = dynVis.getFaceValue(i1);
            //}
        }
        protected void ClearWallPara(int taskIndex)
        {
            foreach (var b in VelVector.NoSlipWalls)
            {
                b.ThreadIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
                for (; i1 < endp1; i1++)
                {
                    Interlocked.Exchange(ref cellCommon[Faces[i1].OSideCell], 0);// = 0;
                }
            }
        }
        public void Add_k_Ab_byConstitutiveEquation(int taskIndex)
        {
            //1.清空数组，准备
            ClearWallPara(taskIndex);
            MultiThreads.Sync(taskIndex);
            //2.计算边界Pk
            {
                foreach (var b in VelVector.NoSlipWalls)
                {
                    b.ThreadIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
                    for (; i1 < endp1; i1++)
                    {
                        ref FvFace face = ref Faces[i1];
                        double Pk = calBound_Pk(i1, ref face);
                        ThreadsMath.KeepMax(ref cellCommon[face.OSideCell], Pk);// * cells[face.OwnerCellIndex].Volumn);
                    }
                }
                MultiThreads.Sync(taskIndex);
            }
            //3.将源项增加的Matrix
            {
                InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
                MatrixForTurbulence mtool = new MatrixForTurbulence();
                for (; i1 < endp1; i1++)
                {
                    double volumn = Cells[i1].Volumn;
                    double des = density.GetCellValue(i1);
                    double w = w_x[i1];
                    //判断是否边界临近
                    //1.如果不紧邻壁面，计算Pk
                    //2.如果紧邻壁面，先设置0，后
                    if (isWallAdj[i1])
                    {
                        k_b[i1] += cellCommon[i1] * volumn;
                    }
                    else
                    {
                        double Pk = mtool.CalPk_Imcompressible(ref U_grad[i1], ref V_grad[i1], ref W_grad[i1], des, k_x[i1], w);
                        //if (i1 == 1221) Console.WriteLine($"Pk nei:{Pk}");
                        k_b[i1] += Pk * volumn;
                    }
                    k_A_Diag[i1] += betastar * des * w * volumn;
                }
            }
        }
        public void Cal_wGamma(int taskIndex)
        {
            cal_kwVis(TotalIndex, sigma, taskIndex);
            MultiThreads.Sync(taskIndex);
            calGamma_byAverage(InnerIndex, taskIndex);
            foreach (var b in VelVector.VelocityInlets)
            {
                calGamma_byBoundNeighbor(b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelVector.PressureInOuts)
            {
                calGamma_byBoundNeighbor(b.ThreadIndex, taskIndex);
            }
            //测试用，先直接选用层流参数
            //totalIndex.getFaceSection(out int i1, out int endp1, taskIndex);
            //for (; i1 < endp1; i1++)
            //{
            //    GammaCommon[i1] = dynVis.getFaceValue(i1);
            //}
        }
        public void Set_w_Ab_byConstitutiveEquation_AfterSetAll(int taskIndex)
        {
            //1.清空数组，准备
            ClearWallPara(taskIndex);
            MultiThreads.Sync(taskIndex);
            //2.计算边界omegaC
            {
                foreach (var b in VelVector.NoSlipWalls)
                {
                    b.ThreadIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
                    for (; i1 < endp1; i1++)
                    {
                        ref FvFace face = ref Faces[i1];
                        double omegaC = calBound_omegaC(i1, ref Faces[i1]);
                        ThreadsMath.KeepMax(ref cellCommon[face.OSideCell], omegaC);
                    }
                }
                MultiThreads.Sync(taskIndex);
            }
            //3.增加源项或指定omegaC
            {
                InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
                MatrixForTurbulence mtool = new MatrixForTurbulence();
                RotationTensor Omegaij = new RotationTensor();
                for (; i1 < endp1; i1++)
                {
                    if (isWallAdj[i1])
                    {
                        double omegaC = cellCommon[i1];
                        w_b[i1] = omegaC * w_A_Diag[i1];
                        w_x[i1] = omegaC;
                        ref FvCell ocell = ref Cells[i1];
                        foreach (int f in ocell.OFaces)
                        {
                            w_A_LU_O[f] = 0;
                        }
                        foreach (int f in ocell.NFaces)
                        {
                            w_A_LU_N[f] = 0;
                        }
                    }
                    else
                    {
                        double volumn = Cells[i1].Volumn;
                        double des = density.GetCellValue(i1);
                        double k = k_x[i1];
                        double w = w_x[i1];
                        ref Vector3d gradU = ref U_grad[i1];
                        ref Vector3d gradV = ref V_grad[i1];
                        ref Vector3d gradW = ref W_grad[i1];
                        double Pk = mtool.CalPk_Imcompressible(ref gradU, ref gradV, ref gradW, des, k, w);
                        double sigmaddkdw = Vector3d.DotProduct(k_grad[i1], w_grad[i1]);
                        if (sigmaddkdw > 0)
                        {
                            sigmaddkdw *= sigmado * des / w;
                        }
                        else
                        {
                            sigmaddkdw = 0;
                        }
                        w_b[i1] += (alpha * w / k * Pk + sigmaddkdw) * volumn;
                        double OOS = Omegaij.cal_OmegaOmegaSijjian(ref gradU, ref gradV, ref gradW, ref mtool);
                        double Chiw = Math.Abs(OOS / Math.Pow(betastar * w, 3));
                        double beta = betao * ((1 + 85 * Chiw) / (1 + 100 * Chiw));
                        w_A_Diag[i1] += beta * des * w * volumn;// rhowvolumn;
                    }
                }
            }
        }

        protected void cal_kwVis(IndexNumber indexNumber, double sigmaStarForkAndsigmaForw, int taskIndex)
        {
            indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                cellCommon[i1] = sigmaStarForkAndsigmaForw * density.GetCellValue(i1) * k_x[i1] / w_x[i1];
            }
        }

        protected void cal_muT_innerCell(int taskIndex)
        {
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            MatrixForTurbulence sij = new MatrixForTurbulence();
            for (; i1 < endp1; i1++)
            {
                cellCommon[i1] = sij.CalTurVis(ref U_grad[i1], ref V_grad[i1], ref W_grad[i1], density.GetCellValue(i1), k_x[i1], w_x[i1]);
                //这个公式待确认，紧邻边界的网格也是这样算吗？！！！！！！！
                //边界处的速度梯度很大，计算出的粘度会很大
            }
        }
        protected void calGamma_byAverage(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                GammaCommon[i1] = dynVis.GetFaceValue(i1) + Faces[i1].CalAverage(cellCommon);
            }
        }
        protected void calGamma_byBoundOwner(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                GammaCommon[i1] = dynVis.GetFaceValue(i1) + cellCommon[Faces[i1].OSideCell];
            }
        }
        protected void calGamma_byBoundNeighbor(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                GammaCommon[i1] = dynVis.GetFaceValue(i1) + cellCommon[Faces[i1].NSideCell];
            }
        }
        protected void calyrhomu(out double y, out double rho, out double mu, out double yrho_mu, int faceIndex, ref FvFace face)
        {
            rho = density.GetCellValue(face.OSideCell);
            mu = dynVis.GetFaceValue(faceIndex);
            y = WallDis_x[face.OSideCell];//第一层cell
            yrho_mu = y * rho / mu;
        }
        protected void calRedgamma(out double gamma, out double kSqrted, ref FvFace face, double yrho_mu)
        {
            kSqrted = Math.Sqrt(k_x[face.OSideCell]);
            double Red = kSqrted * yrho_mu;//这个地方用的Face上的Gamma，为了统一，Gamma都是用的Face
            gamma = Math.Exp(-Red / 11.0);
        }
        protected void calustar(out double uStar, out double u, out double yPlus, ref FvFace face, double gamma, double yrho_mu, double kSqrted)
        {
            Vector3d vel = VelVector.GetVelocity(face.OSideCell) - VelVector.GetVelocity(face.NSideCell);
            vel.MinusToSelf(Vector3d.DotProduct(vel, face.Normal) / face.Area * face.Normal);
            u = vel.CalLength();
            uStar = gamma * Math.Sqrt(u / yrho_mu) + (1 - gamma) * Cmu14 * kSqrted;//感觉应该只包括后一项
            yPlus = yrho_mu * uStar;
        }
        protected double calBound_mu_w(int faceIndex, ref FvFace face)
        {
            calyrhomu(out double y, out double rho, out double mu, out double yrho_mu, faceIndex, ref face);
            calRedgamma(out double gamma, out double kSqrted, ref face, yrho_mu);
            calustar(out double uStar, out double u, out double yPlus, ref face, gamma, yrho_mu, kSqrted);
            //double uPlus = gamma * yPlus + gamma_1 * (Karman_Reciprocal * Math.Log(yPlus) + B);//必须采用更详细的公式，这样弄第二项可能为负//也用了混合函数，StarCCM+帮助中，由更详细的公式，后续可参考改进
            double uPlus = yPlus < 11 ? yPlus : Karman_Reciprocal * Math.Log(yPlus) + B;
            return gamma * mu + (1 - gamma) * (rho * uStar * y) / uPlus;//壁面处粘性系数mu
        }
        protected double calBound_Pk(int faceIndex, ref FvFace face)
        {
            calyrhomu(out double y, out double rho, out double mu, out double yrho_mu, faceIndex, ref face);
            calRedgamma(out double gamma, out double kSqrted, ref face, yrho_mu);
            calustar(out double uStar, out double u, out double yPlus, ref face, gamma, yrho_mu, kSqrted);
            return gamma * mu * Math.Pow(u / y, 2) + (1 - gamma) * rho * Math.Pow(uStar, 3) / (Karman * y);
        }
        protected double calBound_omegaC(int faceIndex, ref FvFace face)
        {
            calyrhomu(out double y, out double rho, out double mu, out double yrho_mu, faceIndex, ref face);
            calRedgamma(out double gamma, out double kSqrted, ref face, yrho_mu);
            return gamma * 60 * mu / (beta1 * rho * Math.Pow(y, 2)) + (1 - gamma) * (kSqrted / (Karman * Cmu14 * y));
        }

        public struct MatrixForTurbulence
        {
            public double a11;
            public double a21;
            public double a22;
            public double a31;
            public double a32;
            public double a33;
            public double CalPk_Imcompressible(ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW, double density, double k, double w)
            {
                double muT = CalTurVis(ref gradU, ref gradV, ref gradW, density, k, w);
                return muT * (2 * (Math.Pow(gradU.X, 2) + Math.Pow(gradV.Y, 2) + Math.Pow(gradW.Z, 2))
                    + Math.Pow(gradU.Y + gradV.X, 2) + Math.Pow(gradU.Z + gradW.X, 2) + Math.Pow(gradV.Z + gradW.Y, 2));
            }
            public double CalPk(ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW, double density, double k, double w)
            {
                double muT = CalTurVis(ref gradU, ref gradV, ref gradW, density, k, w);
                return getPk(ref gradU, ref gradV, ref gradW, muT, density * k);
            }
            public double CalTurVis(ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW, double density, double k, double w)//即muT
            {
                setSij_Sijba_Sijjian_imcompressable(ref gradU, ref gradV, ref gradW);
                double wbo = Math.Max(w, Clim * Math.Sqrt(div_2_betastar * getSijSij()));
                return density * k / wbo;
            }
            private void setSij_Sijba_Sijjian_imcompressable(ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW)//论文公式5，不可压忽略后一项,即等于Sij
            {
                a11 = gradU.X;
                a22 = gradV.Y;
                a33 = gradW.Z;
                a21 = 0.5 * (gradV.X + gradU.Y);
                a31 = 0.5 * (gradU.Z + gradW.X);
                a32 = 0.5 * (gradW.Y + gradV.Z);
            }
            private void setSijba_compressable(ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW)//论文公式5
            {
                setSij_Sijba_Sijjian_imcompressable(ref gradU, ref gradV, ref gradW);
                double part2 = -LxTool.OneThird * (gradU.X + gradV.Y + gradW.Z);
                a11 += part2;
                a22 += part2;
                a33 += part2;
            }
            private void setSijjian_compressable(ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW)//论文公式13
            {
                setSij_Sijba_Sijjian_imcompressable(ref gradU, ref gradV, ref gradW);
                double part2 = -0.5 * (gradU.X + gradV.Y + gradW.Z);
                a11 += part2;
                a22 += part2;
                a33 += part2;
            }
            private double getSijSij()
            {
                return Math.Pow(a11, 2) + Math.Pow(a22, 2) + Math.Pow(a33, 2)
                    + 2 * Math.Pow(a21, 2) + 2 * Math.Pow(a31, 2) + 2 * Math.Pow(a32, 2);
            }
            private double getPk(ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW, double uT, double rhok)
            {
                //论文公式4 计算rho*tij
                // uT *= 2;//变成了2uT
                // rhok *= -2 * MathTool.onethird;//变成了-2/3*rho*k
                //rhok *= -MathTool.onethird;
                //a11 = uT * a11 + rhok;
                //a22 = uT * a22 + rhok;
                //a33 = uT * a33 + rhok;
                //a21 *= uT;
                //a31 *= uT;
                //a32 *= uT;
                ////计算Pk
                //return a11 * gradU.X + a21 * gradU.Y + a31 * gradU.Z
                //    + a21 * gradV.X + a22 * gradV.Y + a32 * gradV.Z
                //    + a31 * gradW.X + a32 * gradW.Y + a33 * gradW.Z;
                //对以上过程简化
                return 2 * uT * (a11 * gradU.X + a21 * gradU.Y + a31 * gradU.Z
                    + a21 * gradV.X + a22 * gradV.Y + a32 * gradV.Z
                    + a31 * gradW.X + a32 * gradW.Y + a33 * gradW.Z)
                    - LxTool.TwoThird * rhok * (gradU.X + gradV.Y + gradW.Z);
            }


        }
        public struct RotationTensor
        {
            public double a21;
            public double a31;
            public double a32;//注意，数据按笔记上的形式组织，只存储下三角矩阵
            //public void setRotationTensor(ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW)
            //{
            //    a21 = 0.5 * (gradV.X - gradU.Y);
            //    a31 = 0.5 * (-gradU.Z + gradW.X);
            //    a32 = 0.5 * (gradW.Y + gradV.Z);
            //}
            public double cal_OmegaOmegaSijjian(ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW, ref MatrixForTurbulence Sijjian)
            {
                a21 = 0.5 * (gradV.X - gradU.Y);
                a31 = 0.5 * (-gradU.Z + gradW.X);
                a32 = 0.5 * (gradW.Y + gradV.Z);
                return -Math.Pow(a21, 2) * (Sijjian.a11 + Sijjian.a22) - Math.Pow(a31, 2) * (Sijjian.a11 + Sijjian.a33)
                   - Math.Pow(a32, 2) * (Sijjian.a33 + Sijjian.a22)
                   + 2 * a21 * a32 * Sijjian.a31 - 2 * a31 * a32 * Sijjian.a21 - 2 * a21 * a31 * Sijjian.a32;
                //Sijjian.setSij_Sijba_Sijjian_imcompressable(ref gradU, ref gradV, ref gradW);
                //double OOS = -Math.Pow(Omegaij.a21, 2) * (Sijjian.a11 + Sijjian.a22) - Math.Pow(Omegaij.a31, 2) * (Sijjian.a11 + Sijjian.a33)
                //    - Math.Pow(Omegaij.a32, 2) * (Sijjian.a33 + Sijjian.a22)
                //    + 2 * Omegaij.a21 * Omegaij.a32 * Sijjian.a31 - 2 * Omegaij.a31 * Omegaij.a32 * Sijjian.a21 - 2 * Omegaij.a21 * Omegaij.a31 * Sijjian.a32;
            }
        }
        public class BoundkValue : Value1Base<double>
        {
            public double TurbIntensitySquared_div_2;
            protected VelocityVector VelVector;
            public void SetParameter(double turbIntensity, VelocityVector VelVector)
            {
                this.TurbIntensitySquared_div_2 = 0.5 * Math.Pow(turbIntensity, 2);
                this.VelVector = VelVector;
            }
            public override double GetCellValue(int cellIndex)
            {
                return TurbIntensitySquared_div_2 * (VelVector.GetVelocity(cellIndex)).CalLengthSquared();
            }
        }
        public class BoundwValue : Value1Base<double>
        {
            double VisRatio;
            double[] k_x;
            FvCell[] cells;
            Value1Base<double> density;
            Value1Base<double> dynVis;
            public void SetParameter(double VisRatio, double[] k_x, FvCell[] cells, Value1Base<double> density, Value1Base<double> dynVis)
            {
                this.VisRatio = VisRatio;
                this.k_x = k_x;
                this.cells = cells;
                this.density = density;
                this.dynVis = dynVis;
            }
            public override double GetCellValue(int cellIndex)
            {
                return density.GetCellValue(cellIndex) * k_x[cellIndex] / dynVis.GetCellValue(cellIndex) * VisRatio;// dynVis.getFaceValue(cells[cellIndex].CellFacesAsN[0])
            }
        }


    }
    public partial class TurbulenceModelkw
    {
        public double PrT;
        public double div_1_PrT;
        public Value1Base<double> ThermalConductivity;
        public Value1Base<double> SpecificHeat;
        public void SetZone_part3_TemperatureGamma(bool isEnableTemperature, Value1Base<double> specificHeat, Value1Base<double> thermalConductivity, double PrT = 8 / 9.0)
        {
            if (!isEnableTemperature)
            {
                this.SpecificHeat = null;
                this.ThermalConductivity = null;
                return;
            }
            this.SpecificHeat = specificHeat;
            this.ThermalConductivity = thermalConductivity;
            this.PrT = PrT;
            this.div_1_PrT = 1 / PrT;
        }
        public void Cal_TemperatureGamma(int taskIndex)
        {
            {
                InnerIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
                for (; i1 < endp1; i1++)
                {
                    GammaCommon[i1] = ThermalConductivity.GetFaceValue(i1) + SpecificHeat.GetFaceValue(i1) * div_1_PrT * Faces[i1].CalAverage(cellCommon);
                }
            }
            foreach (var b in VelVector.NoSlipWalls)
            {
                b.ThreadIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
                for (; i1 < endp1; i1++)
                {
                    GammaCommon[i1] = calBound_ThermalConductivity(i1, ref Faces[i1]);
                }
            }
            foreach (var b in VelVector.Symmetries)
            {
                b.ThreadIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
                for (; i1 < endp1; i1++)
                {
                    GammaCommon[i1] = 0;
                }
            }
            foreach (var b in VelVector.VelocityInlets)
            {
                calBound_ThermalConductivity(b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelVector.PressureInOuts)
            {
                calBound_ThermalConductivity(b.ThreadIndex, taskIndex);
            }
        }
        protected void calBound_ThermalConductivity(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                GammaCommon[i1] = ThermalConductivity.GetFaceValue(i1) + SpecificHeat.GetFaceValue(i1) * div_1_PrT * cellCommon[Faces[i1].OSideCell];
            }
        }
        protected double calBound_ThermalConductivity(int faceIndex, ref FvFace face)
        {
            calyrhomu(out double y, out double rho, out double mu, out double yrho_mu, faceIndex, ref face);
            calRedgamma(out double gamma, out double kSqrted, ref face, yrho_mu);
            calustar(out double uStar, out double u, out double yPlus, ref face, gamma, yrho_mu, kSqrted);
            double specificHeat = SpecificHeat.GetCellValue(face.OSideCell);
            double Pr = specificHeat * mu / ThermalConductivity.GetCellValue(face.OSideCell);
            double Gamma = 0.01 * Math.Pow(Pr * yPlus, 4) / (1 + 5 * Math.Pow(Pr, 3) * yPlus);
            double exp_Gamma = Math.Exp(-Gamma);
            double Pr_div_Prt = Pr * div_1_PrT;
            double P = 9.24 * (Math.Pow(Pr_div_Prt, 0.75) - 1) * (1 + 0.28 * Math.Exp(-0.007 * Pr_div_Prt));
            double TPlus = exp_Gamma * Pr * yPlus + (1 - exp_Gamma) * PrT * (Karman_Reciprocal * Math.Log(yPlus) + P);
            return rho * specificHeat * uStar * y / TPlus;
        }

    }
}
//protected double calBound_muw(int faceidx)
//{
//    ref Mesh.Face face = ref faces[faceidx];
//    int ocellidx = face.OwnerCellIndex;//第一层cell
//    double rho = density.getCellValue(ocellidx);
//    double mu = dynVis.getFaceValue(faceidx);
//    double y = WallDis_x[ocellidx];
//    double yrho_mu = y * rho / mu;

//    double kSqrted = Math.Sqrt(k_x[ocellidx]);
//    double Red = kSqrted * yrho_mu;//这个地方用的Face上的Gamma，为了统一，Gamma都是用的Face
//    double gamma = Math.Exp(-Red / 11.0);
//    double gamma_1 = 1 - gamma;

//    Vector3d vel = VelVector.getVelocity(ocellidx) - VelVector.getVelocity(face.NeighborCellIndex);
//    vel = vel - Vector3d.DotProduct(vel, face.Normal) / face.Area * face.Normal;
//    double u = vel.CalLength();
//    double uStar = gamma * Math.Sqrt(u / yrho_mu) + gamma_1 * Cmu14 * kSqrted;//感觉应该只包括后一项
//    double yPlus = yrho_mu * uStar;

//    //double uPlus = gamma * yPlus + gamma_1 * (Karman_Reciprocal * Math.Log(yPlus) + B);//必须采用更详细的公式，这样弄第二项可能为负//也用了混合函数，StarCCM+帮助中，由更详细的公式，后续可参考改进
//    double uPlus = yPlus < 11 ? yPlus : Karman_Reciprocal * Math.Log(yPlus) + B;
//    //double muw =
//    return gamma * mu + gamma_1 * (rho * uStar * y) / uPlus;//壁面处粘性系数mu

//    //double Pk = gamma * mu * Math.Pow(u / y, 2) + gamma_1 * rho * Math.Pow(uStar, 3) / (Karman * y);
//    //double omegaC = gamma * 6 * mu / (beta1 * rho * Math.Pow(y, 2)) + gamma_1 * (kSqrted / (Karman * Cmu14 * y));
//}
//protected double calBound_Pk(int faceidx, ref Mesh.Face face)
//{
//    int ocellidx = face.OwnerCellIndex;//第一层cell
//    double rho = density.getCellValue(ocellidx);
//    double mu = dynVis.getFaceValue(faceidx);
//    double y = WallDis_x[ocellidx];
//    double yrho_mu = y * rho / mu;

//    double kSqrted = Math.Sqrt(k_x[ocellidx]);
//    double Red = kSqrted * yrho_mu;//这个地方用的Face上的Gamma，为了统一，Gamma都是用的Face
//    double gamma = Math.Exp(-Red / 11.0);
//    double gamma_1 = 1 - gamma;

//    Vector3d vel = VelVector.getVelocity(ocellidx) - VelVector.getVelocity(face.NeighborCellIndex);
//    vel = vel - Vector3d.DotProduct(vel, face.Normal) / face.Area * face.Normal;
//    double u = vel.CalLength();
//    double uStar = gamma * Math.Sqrt(u / yrho_mu) + gamma_1 * Cmu14 * kSqrted;//感觉应该只包括后一项
//                                                                              //double yPlus = yrho_mu * uStar;
//                                                                              //double uPlus = gamma * yPlus + gamma_1 * (Karman_Reciprocal * Math.Log(yPlus) + B);//也用了混合函数，StarCCM+帮助中，由更详细的公式，后续可参考改进

//    //double muw = gamma * mu + gamma_1 * (rho * uStar * y) / uPlus;//壁面处粘性系数mu

//    //double Pk = gamma * mu * Math.Pow(u / y, 2) + gamma_1 * rho * Math.Pow(uStar, 3) / (Karman * y);
//    return gamma * mu * Math.Pow(u / y, 2) + gamma_1 * rho * Math.Pow(uStar, 3) / (Karman * y);
//    //double omegaC = gamma * 6 * mu / (beta1 * rho * Math.Pow(y, 2)) + gamma_1 * (kSqrted / (Karman * Cmu14 * y));
//}
//protected double calBound_omegaC(int faceidx)
//{
//    ref Mesh.Face face = ref faces[faceidx];
//    int ocellidx = face.OwnerCellIndex;//第一层cell
//    double rho = density.getCellValue(ocellidx);
//    double mu = dynVis.getFaceValue(faceidx);
//    double y = WallDis_x[ocellidx];
//    double yrho_mu = y * rho / mu;

//    double kSqrted = Math.Sqrt(k_x[ocellidx]);
//    double Red = kSqrted * yrho_mu;//这个地方用的Face上的Gamma，为了统一，Gamma都是用的Face
//    double gamma = Math.Exp(-Red / 11.0);
//    double gamma_1 = 1 - gamma;

//    //Vector3d vel = VelVector.getVelocity(ocellidx) - VelVector.getVelocity(face.NeighborCellIndex);
//    //vel = vel - Vector3d.DotProduct(vel, face.Normal) / face.Area * face.Normal;
//    //double u = vel.CalLength();
//    //double uStar = gamma * Math.Sqrt(u / yrho_mu) + gamma_1 * Cmu14 * kSqrted;//感觉应该只包括后一项
//    //double yPlus = yrho_mu * uStar;
//    //double uPlus = gamma * yPlus + gamma_1 * (Karman_Reciprocal * Math.Log(yPlus) + B);//也用了混合函数，StarCCM+帮助中，由更详细的公式，后续可参考改进

//    //double muw = gamma * mu + gamma_1 * (rho * uStar * y) / uPlus;//壁面处粘性系数mu

//    //double Pk = gamma * mu * Math.Pow(u / y, 2) + gamma_1 * rho * Math.Pow(uStar, 3) / (Karman * y);
//    //double omegaC = gamma * 6 * mu / (beta1 * rho * Math.Pow(y, 2)) + gamma_1 * (kSqrted / (Karman * Cmu14 * y));
//    return gamma * 60 * mu / (beta1 * rho * Math.Pow(y, 2)) + gamma_1 * (kSqrted / (Karman * Cmu14 * y));
//}
//public struct BoundParaCalulator
//{
//    Value1Base density;
//    Value1Base dynVis;
//    double[] WallDis_x;
//    double[] k_x;
//    VelocityVector VelVector;
//    public void SetCommon(Value1Base density, Value1Base dynVis, double[] WallDis_x, double[] k_x, VelocityVector VelVector)
//    {
//        this.density = density;
//        this.dynVis = dynVis;
//        this.WallDis_x = WallDis_x;
//        this.k_x = k_x;
//        this.VelVector = VelVector;
//    }
//    int ocellidx;
//    double rho;
//    double mu;
//    double y;
//    double yrho_mu;
//    private void cal_yrhomu(int faceidx, ref Mesh.Face face)
//    {
//        ocellidx = face.OwnerCellIndex;//第一层cell
//        rho = density.getCellValue(ocellidx);
//        mu = dynVis.getFaceValue(faceidx);
//        y = WallDis_x[ocellidx];
//        yrho_mu = y * rho / mu;
//    }
//    double kSqrted;
//    double Red;
//    double gamma;
//    double gamma_1;
//    private void calgamma()
//    {
//        kSqrted = Math.Sqrt(k_x[ocellidx]);
//        Red = kSqrted * yrho_mu;//这个地方用的Face上的Gamma，为了统一，Gamma都是用的Face
//        gamma = Math.Exp(-Red / 11.0);
//        gamma_1 = 1 - gamma;
//    }
//    Vector3d vel;
//    double u;
//    double uStar;
//    private void calustar(ref Mesh.Face face)
//    {
//        vel = VelVector.getVelocity(ocellidx) - VelVector.getVelocity(face.NeighborCellIndex);
//        vel = vel - Vector3d.DotProduct(vel, face.Normal) / face.Area * face.Normal;
//        u = vel.CalLength();
//        uStar = gamma * Math.Sqrt(u / yrho_mu) + gamma_1 * Cmu14 * kSqrted;//感觉应该只包括后一项
//    }
//    double yPlus;
//    private void calyplus()
//    {
//        yPlus = yrho_mu * uStar;
//    }
//    public double Cal_mu_w(int faceidx, ref Mesh.Face face)
//    {
//        cal_yrhomu(faceidx, ref face);
//        calgamma();
//        calustar(ref face);
//        calyplus();
//        //double uPlus = gamma * yPlus + gamma_1 * (Karman_Reciprocal * Math.Log(yPlus) + B);//必须采用更详细的公式，这样弄第二项可能为负//也用了混合函数，StarCCM+帮助中，由更详细的公式，后续可参考改进
//        double uPlus = yPlus < 11 ? yPlus : Karman_Reciprocal * Math.Log(yPlus) + B;
//        return gamma * mu + gamma_1 * (rho * uStar * y) / uPlus;//壁面处粘性系数mu
//    }
//    public double Cal_Pk(int faceidx, ref Mesh.Face face)
//    {
//        cal_yrhomu(faceidx, ref face);
//        calgamma();
//        calustar(ref face);
//        return gamma * mu * Math.Pow(u / y, 2) + gamma_1 * rho * Math.Pow(uStar, 3) / (Karman * y);
//    }
//    public double Cal_omegaC(int faceidx, ref Mesh.Face face)
//    {
//        cal_yrhomu(faceidx, ref face);
//        calgamma();
//        return gamma * 60 * mu / (beta1 * rho * Math.Pow(y, 2)) + gamma_1 * (kSqrted / (Karman * Cmu14 * y));
//    }
//    public double Cal_ThermalConductivity(int faceidx, ref Mesh.Face face, double Pr, double Prt, double specificHeat)
//    {
//        cal_yrhomu(faceidx, ref face);
//        calgamma();
//        calustar(ref face);
//        calyplus();
//        double Gamma = 0.01 * Math.Pow(Pr * yPlus, 4) / (1 + 5 * Math.Pow(Pr, 3) * yPlus);
//        double exp_Gamma = Math.Exp(-Gamma);
//        double Pr_div_Prt = Pr / Prt;
//        double P = 9.24 * (Math.Pow(Pr_div_Prt, 0.75) - 1) * (1 + 0.28 * Math.Exp(-0.007 * Pr_div_Prt));
//        double TPlus = exp_Gamma * Pr * yPlus + (1 - exp_Gamma) * Prt * (Karman_Reciprocal * Math.Log(yPlus) + P);
//        return rho * specificHeat * uStar * y / TPlus;
//    }
//}
//public class ValuePrLReciprocal : Value1Base
//{
//    Value1Base specificHeat;
//    Value1Base dynamicViscosity;
//    Value1Base thermalConductivity;
//    public ValuePrLReciprocal(Value1Base specificHeat, Value1Base dynamicViscosity, Value1Base thermalConductivity)
//    {
//        this.specificHeat = specificHeat;
//        this.dynamicViscosity = dynamicViscosity;
//        this.thermalConductivity = thermalConductivity;
//    }
//    public override double getCellValue(int cellIndex)
//    {
//        return thermalConductivity.getCellValue(cellIndex) / (specificHeat.getCellValue(cellIndex) * dynamicViscosity.getCellValue(cellIndex));
//    }
//    public override double getFaceValue(int faceIndex)
//    {
//        return thermalConductivity.getFaceValue(faceIndex) / (specificHeat.getFaceValue(faceIndex) * dynamicViscosity.getFaceValue(faceIndex));
//    }
//}
//public TurbulentViscosityOfkw TauForMom;
//protected double[] TurVisOnFace;//还有k w方程的//湍流粘度，muT////rho*k/w~论文公式6
//public TurbulentViscosityOfkw TauForkw;//k w方程公用一个Tau，注意切换，先算w后算k，能量方程的系数好像与k方程的相同
//protected double[] TauForkw_xf;
//protected double[] rhokw;//rho*k/w

//totalIndex.getCellSection(out int i1, out int endp1, taskIndex);
//for (; i1 < endp1; i1++)
//{
//    rhokw[i1] = sigma * density.getCellValue(i1) * k_x[i1] / w_x[i1];
//}
//protected void updateTauForkw_innerCell(int taskIndex)
//{
//    innerIndex.getCellSection(out int i1, out int endp1, taskIndex);
//    for (; i1 < endp1; i1++)
//    {
//        rhokw[i1] = density.getCellValue(i1) * k_x[i1] / w_x[i1];
//    }
//}
//protected void updateTauFork_innerFace(int taskIndex)
//{
//    updateTauForkw_innerFace(sigmastar, taskIndex);
//}
//protected void updateTauForw_innerFace(int taskIndex)
//{
//    updateTauForkw_innerFace(sigma, taskIndex);
//}
//protected void updateTauForkw_innerFace(double a, int taskIndex)
//{
//    innerIndex.getFaceSection(out int i1, out int endp1, taskIndex);
//    for (; i1 < endp1; i1++)
//    {
//        TauForkw_xf[i1] = dynVis.getFaceValue(i1) + a * faces[i1].calAverageValue(rhokw);
//    }
//}
//public struct SymMatrix3
//{
//    public double a11;
//    public double a21;
//    public double a22;
//    public double a31;
//    public double a32;
//    public double a33;
//    public void setSij_Sijba_Sijjian_imcompressable(ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW)//论文公式5，不可压忽略后一项,即等于Sij
//    {
//        a11 = gradU.X;
//        a22 = gradV.Y;
//        a33 = gradW.Z;
//        a21 = 0.5 * (gradV.X + gradU.Y);
//        a31 = 0.5 * (gradU.Z + gradW.X);
//        a32 = 0.5 * (gradW.Y + gradV.Z);
//    }
//    public void setSijba_compressable(ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW)//论文公式5
//    {
//        setSij_Sijba_Sijjian_imcompressable(ref gradU, ref gradV, ref gradW);
//        double part2 = -MathTool.onethird * (gradU.X + gradV.Y + gradW.Z);
//        a11 += part2;
//        a22 += part2;
//        a33 += part2;
//    }
//    public void setSijjian_compressable(ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW)//论文公式13
//    {
//        setSij_Sijba_Sijjian_imcompressable(ref gradU, ref gradV, ref gradW);
//        double part2 = -0.5 * (gradU.X + gradV.Y + gradW.Z);
//        a11 += part2;
//        a22 += part2;
//        a33 += part2;
//    }
//    public double getSijSij()
//    {
//        return Math.Pow(a11, 2) + Math.Pow(a22, 2) + Math.Pow(a33, 2)
//            + 2 * Math.Pow(a21, 2) + 2 * Math.Pow(a31, 2) + 2 * Math.Pow(a32, 2);
//    }

//    public double calPk_imcompressible(ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW, double uT, double rhok)
//    {
//        setSij_Sijba_Sijjian_imcompressable(ref gradU, ref gradV, ref gradW);
//        //论文公式4 计算rho*tij
//        // uT *= 2;//变成了2uT
//        // rhok *= -2 * MathTool.onethird;//变成了-2/3*rho*k
//        //rhok *= -MathTool.onethird;
//        //a11 = uT * a11 + rhok;
//        //a22 = uT * a22 + rhok;
//        //a33 = uT * a33 + rhok;
//        //a21 *= uT;
//        //a31 *= uT;
//        //a32 *= uT;
//        ////计算Pk
//        //return a11 * gradU.X + a21 * gradU.Y + a31 * gradU.Z
//        //    + a21 * gradV.X + a22 * gradV.Y + a32 * gradV.Z
//        //    + a31 * gradW.X + a32 * gradW.Y + a33 * gradW.Z;
//        //对以上过程简化
//        return 2 * uT * (a11 * gradU.X + a21 * gradU.Y + a31 * gradU.Z
//            + a21 * gradV.X + a22 * gradV.Y + a32 * gradV.Z
//            + a31 * gradW.X + a32 * gradW.Y + a33 * gradW.Z)
//            - MathTool.twothird * rhok * (gradU.X + gradV.Y + gradW.Z);
//    }
//    public double getPk(ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW, double uT, double rhok)
//    {
//        return 2 * uT * (a11 * gradU.X + a21 * gradU.Y + a31 * gradU.Z
//            + a21 * gradV.X + a22 * gradV.Y + a32 * gradV.Z
//            + a31 * gradW.X + a32 * gradW.Y + a33 * gradW.Z)
//            - MathTool.twothird * rhok * (gradU.X + gradV.Y + gradW.Z);
//    }
//}
//protected void calTotalVis_Mom_innerCell(int taskIndex)
//{
//    innerIndex.getCellSection(out int i1, out int endp1, taskIndex);
//    SymMatrix3 sij = new SymMatrix3();
//    double sijxishu = 2 / betastar;
//    for (; i1 < endp1; i1++)
//    {
//        sij.setSij_Sijba_Sijjian_imcompressable(ref U_grad[i1], ref V_grad[i1], ref W_grad[i1]);
//        double wbo = Math.Max(w_x[i1], Clim * Math.Sqrt(sijxishu * sij.getSijSij()));
//        TurVis[i1] = density.getCellValue(i1) * k_x[i1] / wbo;
//    }
//}

//public class TurbulentViscosityOfkw : Value1Base//Taut
//{
//    public TurbulentViscosityOfkw()
//    {

//    }
//    protected double[] x_common;
//    public double[] x_OnFace;
//    public override double getFaceValue(int faceIndex)
//    {
//        return x_OnFace[faceIndex];
//    }
//}
//public void Add_kw_Ab(int taskIndex)
//{
//    innerIndex.getCellSection(out int i1, out int endp1, taskIndex);
//    SymMatrix3 symmatrix = new SymMatrix3();
//    SymMatrix3 Sijjian = new SymMatrix3();
//    RotationTensor Omegaij = new RotationTensor();
//    for (; i1 < endp1; i1++)
//    {
//        double volumn = cells[i1].Volumn;
//        double des = density.getCellValue(i1);
//        double k = k_x[i1];
//        double w = w_x[i1];
//        double rhowvolumn = des * w * volumn;//两处用到
//        ref Vector3d gradU = ref U_grad[i1];
//        ref Vector3d gradV = ref V_grad[i1];
//        ref Vector3d gradW = ref W_grad[i1];
//        double Pk = symmatrix.calPk_imcompressible(ref gradU, ref gradV, ref gradW, TurVis[i1], des * k);
//        k_b[i1] += Pk * volumn;
//        k_A_Diag[i1] += betastar * rhowvolumn;
//        double sigmaddkdw = Vector3d.DotProduct(k_grad[i1], w_grad[i1]);
//        if (sigmaddkdw > 0)
//        {
//            sigmaddkdw *= sigmado * des / w;
//        }
//        else
//        {
//            sigmaddkdw = 0;
//        }
//        w_b[i1] += (alpha * w / k * Pk + sigmaddkdw) * volumn;
//        Omegaij.setRotationTensor(ref gradU, ref gradV, ref gradW);
//        Sijjian.setSij_Sijba_Sijjian_imcompressable(ref gradU, ref gradV, ref gradW);
//        double OOS = -Math.Pow(Omegaij.a21, 2) * (Sijjian.a11 + Sijjian.a22) - Math.Pow(Omegaij.a31, 2) * (Sijjian.a11 + Sijjian.a33)
//            - Math.Pow(Omegaij.a32, 2) * (Sijjian.a33 + Sijjian.a22)
//            + 2 * Omegaij.a21 * Omegaij.a32 * Sijjian.a31 - 2 * Omegaij.a31 * Omegaij.a32 * Sijjian.a21 - 2 * Omegaij.a21 * Omegaij.a31 * Sijjian.a32;
//        double Chiw = Math.Abs(OOS / Math.Pow(betastar * w, 3));
//        double beta = betao * ((1 + 85 * Chiw) / (1 + 100 * Chiw));
//        w_A_Diag[i1] += beta * rhowvolumn;
//    }
//}
//protected double calChiw(ref RotationTensor Omegaij, ref SymMatrix3 Sijjian, double w)
//{
//    double fenzi = -Math.Pow(Omegaij.a21, 2) * (Sijjian.a11 + Sijjian.a22)
//        - Math.Pow(Omegaij.a31, 2) * (Sijjian.a11 + Sijjian.a33)
//        - Math.Pow(Omegaij.a32, 2) * (Sijjian.a33 + Sijjian.a22)
//        + 2 * Omegaij.a21 * Omegaij.a32 * Sijjian.a31 - 2 * Omegaij.a31 * Omegaij.a32 * Sijjian.a21 - 2 * Omegaij.a21 * Omegaij.a31 * Sijjian.a32;
//    return Math.Abs(fenzi / Math.Pow(betastar * w, 3));
//}
//public void Add_UVWbound_Ab(IndexNumber indexNumber, int taskIndex)
//{
//    indexNumber.getFaceSection(out int faceidx, out int endp1, taskIndex);
//    for (; faceidx < endp1; faceidx++)
//    {
//        ref Mesh.Face face = ref faces[faceidx];
//        ref Mesh.Cell ocell = ref cells[face.OwnerCellIndex];
//        int ocellidx = face.OwnerCellIndex;//第一层cell
//        double rho = density.getCellValue(ocellidx);
//        double mu = dynVis.getFaceValue(faceidx);
//        double y = WDis_x[ocellidx];
//        double yrho_mu = y * rho / mu;

//        double kSqrted = Math.Sqrt(k_x[ocellidx]);
//        double Red = kSqrted * yrho_mu;//这个地方用的Face上的Gamma，为了统一，Gamma都是用的Face
//        double gamma = Math.Exp(-Red / 11.0);
//        double gamma_1 = 1 - gamma;

//        Vector3d vel = VelVector.getVelocity(ocellidx) - VelVector.getVelocity(face.NeighborCellIndex);
//        vel = vel - Vector3d.DotProduct(vel, face.Normal) / face.Area * face.Normal;
//        double u = vel.CalLengthSquared();
//        double uStar = gamma * Math.Sqrt(u / yrho_mu) + gamma_1 * Cmu14 * kSqrted;//感觉应该只包括后一项
//        double yPlus = yrho_mu * uStar;
//        double uPlus = gamma * yPlus + gamma_1 * (Karman_Reciprocal * Math.Log(yPlus) + B);//也用了混合函数，StarCCM+帮助中，由更详细的公式，后续可参考改进

//        double muw = gamma * mu + gamma_1 * (rho * uStar * y) / uPlus;//壁面处粘性系数mu
//        TotalVisMom[faceidx] = muw;

//        double Pk = gamma * mu * Math.Pow(u / y, 2) + gamma_1 * rho * Math.Pow(uStar, 3) / (Karman * y);
//        double omegaC = gamma * 6 * mu / (beta1 * rho * Math.Pow(y, 2)) + gamma_1 * (kSqrted / (Karman * Cmu14 * y));
//    }
//}


//protected double calBound_muw(ref Mesh.Face face, double rho, double mu, double y, double k, Vector3d vel)//这个地方mu用的Face上的Gamma，为了统一，Gamma都是用的Face
//{
//    double yrho_mu = y * rho / mu;

//    double kSqrted = Math.Sqrt(k);
//    double Red = kSqrted * yrho_mu;
//    double gamma = Math.Exp(-Red / 11.0);
//    double gamma_1 = 1 - gamma;

//    vel.MinusToSelf(Vector3d.DotProduct(vel, face.Normal) / face.Area * face.Normal);
//    double u = vel.CalLengthSquared();
//    double uStar = gamma * Math.Sqrt(u / yrho_mu) + gamma_1 * Cmu14 * kSqrted;//感觉应该只包括后一项
//    double yPlus = yrho_mu * uStar;
//    double uPlus = gamma * yPlus + gamma_1 * (Karman_Reciprocal * Math.Log(yPlus) + B);//也用了混合函数，StarCCM+帮助中，由更详细的公式，后续可参考改进

//    double muw = gamma * mu + gamma_1 * (rho * uStar * y) / uPlus;//壁面处粘性系数mu
//    return muw;
//    //double Pk = gamma * mu * Math.Pow(u / y, 2) + gamma_1 * rho * Math.Pow(uStar, 3) / (Karman * y);
//    //double omegaC = gamma * 6 * mu / (beta1 * rho * Math.Pow(y, 2)) + gamma_1 * (kSqrted / (Karman * Cmu14 * y));
//}
//protected double calBound_Pk(ref Mesh.Face face, double rho, double mu, double y, double k, Vector3d vel)//这个地方mu用的Face上的Gamma，为了统一，Gamma都是用的Face
//{
//    double yrho_mu = y * rho / mu;

//    double kSqrted = Math.Sqrt(k);
//    double Red = kSqrted * yrho_mu;
//    double gamma = Math.Exp(-Red / 11.0);
//    double gamma_1 = 1 - gamma;

//    vel.MinusToSelf(Vector3d.DotProduct(vel, face.Normal) / face.Area * face.Normal);
//    double u = vel.CalLengthSquared();
//    double uStar = gamma * Math.Sqrt(u / yrho_mu) + gamma_1 * Cmu14 * kSqrted;
//    //感觉应该只包括后一项
//    //double yPlus = yrho_mu * uStar;
//    // double uPlus = gamma * yPlus + gamma_1 * (Karman_Reciprocal * Math.Log(yPlus) + B);//也用了混合函数，StarCCM+帮助中，由更详细的公式，后续可参考改进

//    // double muw = gamma * mu + gamma_1 * (rho * uStar * y) / uPlus;//壁面处粘性系数mu

//    double Pk = gamma * mu * Math.Pow(u / y, 2) + gamma_1 * rho * Math.Pow(uStar, 3) / (Karman * y);
//    return Pk;
//    // double omegaC = gamma * 6 * mu / (beta1 * rho * Math.Pow(y, 2)) + gamma_1 * (kSqrted / (Karman * Cmu14 * y));
//}
//protected double calBound_omegaC(double rho, double mu, double y, double k)//这个地方mu用的Face上的Gamma，为了统一，Gamma都是用的Face
//{
//    double yrho_mu = y * rho / mu;

//    double kSqrted = Math.Sqrt(k);
//    double Red = kSqrted * yrho_mu;
//    double gamma = Math.Exp(-Red / 11.0);
//    // double gamma_1 = 1 - gamma;

//    // vel.MinusToSelf(Vector3d.DotProduct(vel, face.Normal) / face.Area * face.Normal);
//    // double u = vel.CalLengthSquared();
//    //  double uStar = gamma * Math.Sqrt(u / yrho_mu) + gamma_1 * Cmu14 * kSqrted;//感觉应该只包括后一项
//    //double yPlus = yrho_mu * uStar;
//    // double uPlus = gamma * yPlus + gamma_1 * (Karman_Reciprocal * Math.Log(yPlus) + B);//也用了混合函数，StarCCM+帮助中，由更详细的公式，后续可参考改进

//    // double muw = gamma * mu + gamma_1 * (rho * uStar * y) / uPlus;//壁面处粘性系数mu

//    //double Pk = gamma * mu * Math.Pow(u / y, 2) + gamma_1 * rho * Math.Pow(uStar, 3) / (Karman * y);
//    double omegaC = gamma * 6 * mu / (beta1 * rho * Math.Pow(y, 2)) + (1 - gamma) * (kSqrted / (Karman * Cmu14 * y));
//    return omegaC;
//}
//protected double calTurVis(ref SymMatrix3 MatrixTool, ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW, double density, double k, double w)
//{
//    MatrixTool.setSij_Sijba_Sijjian_imcompressable(ref gradU, ref gradV, ref gradW);
//    double wbo = Math.Max(w, Clim * Math.Sqrt(div_2_betastar * MatrixTool.getSijSij()));
//    return density * k / wbo;
//}
//protected double calPk(ref SymMatrix3 MatrixTool, ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW, double density, double k, double w)
//{
//    double muT = calTurVis(ref MatrixTool, ref gradU, ref gradV, ref gradW, density, k, w);
//    return MatrixTool.calPk_imcompressible(ref gradU, ref gradV, ref gradW, muT, density * k);
//}